program test_solve

    use gauss
    use exp, only: new_solve => solve, solve_f => solve1
    integer :: ipiv(2)

    real(kind=8) :: A(2,2), b(2, 1)
    
    A = reshape([1, 2, 0, 1], [2, 2])
    b = reshape([3, 1      ], [2, 1])

    print *, solve_f(A, b)

    call dgesv( 2, 1, A, 2, ipiv, b, 2, info)

    print *, b

end program test_solve